\section{概述}

在用数学模型, 包括概率统计模型处理实际应用中的问题时, 我们希望建立的模型能够尽可能地符合实际情况。但是, 实际情况是错综复杂的, 如果一味地要求模型与实际完全相符, 会导致模型过于复杂, 以至于不能进行严格理论分析, 结果导致模型不能使用。所以, 实际建模时会忽略许多细节, 增加一些可能很难验证的理论假设, 使得模型比较简单, 可以用数学理论进行分析研究。

这样, 简化的模型就可以与实际情况有较大的差距, 即使我们对模型进行了完美的理论分析, 也不能保证分析结果是可信的。这一困难可以用随机模拟的方法解决。

{\color{red}\textbf{[关键理解]:} 理论分析像``纸上推演棋局''——需要简化规则才能算出最优策略,但简化后的结果未必符合真实棋局。随机模拟则是``实际下棋千万盘''——虽然单盘有随机性,但大量重复能揭示真实规律,绕过理论分析的困难。}

模拟是指把某一现实的或抽象的系统的某种特征或部分状态, 用另一系统 (称为模拟模型)来代替或模拟。为了解决某问题, 把它变成一个概率模型的求解问题, 然后产生符合模型的大量随机数, 对产生的随机数进行分析从而求解问题, 这种方法叫做随机模拟方法, 又称为蒙特卡洛 (Monte Carlo) 方法。

{\color{red}\textbf{[历史背景]:} ``蒙特卡洛''名称来自摩纳哥赌场,由冯·诺依曼和乌拉姆在曼哈顿计划(1940年代)中使用这个代号。核心思想就像赌徒``试手气千万次''——单次输赢随机,但大数定律保证长期统计规律可预测。}

例如, 一个交通路口需要找到一种最优的控制红绿灯信号的办法, 使得通过路口的汽车耽搁的平均时间最短,而行人等候过路的时间不超过某一给定的心理极限值。十字路口的信号共有四个方向, 每个方向又分直行、左转、右转。因为汽车和行人的到来是随机的, 我们要用随机过程来描述四个方向的汽车到来和路口的行人到来过程。理论建模分析很难解决这个最优化问题。但是, 我们可以采集汽车和行人到来的频率, 用随机模拟方法模拟汽车和行人到来的过程, 并模拟各种控制方案, 记录不同方案造成的等待时间, 通过统计比较找出最优的控制方案。

随机模拟中的随机性可能来自模型本身的随机变量, 比如上面描述的汽车和行人到来, 也可能是把非随机的问题转换为概率模型的特征量估计问题从而用随机模拟方法解决。

例 3.1.1. 为了计算圆周率 \(\pi\) 的近似值可以使用随机模拟方法。如果向正方形 \(D = \{ \left( {x,y}\right)\) : \(x \in  \left\lbrack  {-1,1}\right\rbrack  ,y \in  \left\lbrack  {-1,1}\right\rbrack\) 内随机等可能投点,落入单位圆 \(C = \left\{  {\left( {x,y}\right)  : {x}^{2} + {y}^{2} \leq  1}\right\}\) 的概率为面积之比 \(p = \frac{\pi }{4}\) 。如果独立重复地投了 \(N\) 个点,落入 \(C\) 中的点的个数 \(\xi\) 的平均值为 \({E\xi } = {pN}\) , 由概率的频率解释,

\[
\frac{\xi }{N} \approx  \frac{\pi }{4},\;\pi  \approx  \widehat{\pi } = \frac{4\xi }{N}
\]

可以这样给出 \(\pi\) 的近似值。

{\color{red}\textbf{[物理图像]:} 想象向靶子扔飞镖——正方形靶子内嵌一个圆形靶心。闭着眼睛随机扔,命中靶心的比例就是圆与正方形的面积比$\pi/4$。扔得越多,这个比例越接近真实值。这把几何问题(计算$\pi$)转化为统计问题(数飞镖)。}

随机模拟方法会引入所谓随机模拟误差: 上例中估计的 \(\widehat{\pi }\) 实际是随机的,如果再次独立重复投 \(N\) 个点,得到的 \(\widehat{\pi }\) 和上一次结果会有不同。这是随机模拟方法的特点,即结果具有随机性。因为结果的随机性导致的误差叫做随机模拟误差。

{\color{red}\textbf{[解释]:} 这不是``算错了'',而是方法本质——就像测量同一个人的身高10次,每次结果略有不同(站姿、仪器读数误差)。随机模拟的``测量仪器''就是随机数,每次运行都会产生不同的随机序列,因此结果必然有波动。}

使用随机模拟方法, 我们必须了解随机模拟误差的大小, 这样我们才能够设计合适的重复试验次数来控制随机模拟误差。比如,这个例子中 \(\xi\) 服从 \(\mathrm{B}\left( {N,\frac{\pi }{4}}\right)\) 分布,有

\[
\operatorname{Var}\left( \widehat{\pi }\right)  = \frac{\pi \left( {4 - \pi }\right) }{16N},
\]

由中心极限定理, \(\widehat{\pi }\) 近似服从 \(\mathrm{N}\left( {\pi ,\frac{\pi \left( {4 - \pi }\right) }{16N}}\right)\) 分布,所以随机模拟误差的幅度大约在 \(\pm  2\sqrt{\frac{\pi \left( {4 - \pi }\right) }{16N}}\) (随机模拟误差 95\% 以上落入此区间)。

{\color{red}\textbf{[关键理解]:} 注意方差中的$1/N$——标准差是$1/\sqrt{N}$。这意味着要把误差减半,需要4倍样本量;要精度提高10倍,需要100倍计算量。这是所有蒙特卡洛方法的``宿命''——收敛慢但稳定,与问题维度无关。}

这样的随机误差幅度也是随机模拟误差的典型情况。一般地,假设随机变量 \(X\) 期望为 \(\theta\) ,方差为 \({\sigma }^{2}\) 。产生随机变量 \(X\) 的 \(N\) 个独立同分布随机数 \({X}_{i},i = 1,2,\ldots ,N\) ,用样本平均值 \({\bar{X}}_{N} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}{X}_{i}\) 估计 \(\theta\) ,由中心极限定理可知 \({\bar{X}}_{N}\) 近似服从 \(\mathrm{N}\left( {\theta ,{\sigma }^{2}/N}\right)\) ,于是随机模拟误差幅度为 \({O}_{p}\left( \frac{\sigma }{\sqrt{N}}\right) \left( {N \rightarrow  \infty }\right)\) 。为了估计 \({\bar{X}}_{N}\) 的渐近标准差 \(\sigma /\sqrt{N}\) ,可以用样本方差 \({S}_{N}^{2} = \frac{1}{N - 1}\mathop{\sum }\limits_{{i = 1}}^{N}{\left( {X}_{i} - {\bar{X}}_{N}\right) }^{2}\) 代替 \({\sigma }^{2}\) 进行计算。为了控制估计的标准差小于 \({\sigma }_{0}\) ,可以先取较小的 \({N}_{0}\) 抽取 \({N}_{0}\) 个样本值计算出 \({S}_{{N}_{0}}^{2}\) ,用 \({S}_{{N}_{0}}^{2}\) 估计 \({\sigma }^{2}\) ,然后求需要的 \(N\) 的大小:

\[
\frac{{S}_{{N}_{0}}}{\sqrt{N}} < {\sigma }_{0},\;N > \frac{{S}_{{N}_{0}}^{2}}{{\sigma }_{0}^{2}}.
\]

用 \({\bar{X}}_{N}\) 估计 \(\theta  = {EX}\) 时,也可以利用中心极限定理计算 \(\theta\) 的近似 \({95}\%\) 置信区间:

\[
{\bar{X}}_{N} \pm  2{S}_{N}/\sqrt{N}. \tag{3.1}
\]

{\color{red}\textbf{[为什么是$\sqrt{N}$？]:} 这来自独立随机变量求和的根本性质。$N$个独立同分布随机变量的和,方差是单个变量方差的$N$倍,因此标准差是$\sqrt{N}$倍。平均后除以$N$,标准差变为$1/\sqrt{N}$。这是中心极限定理的核心——涨落被``部分''抵消,但不是完全抵消。}

随机模拟方法虽然避免了复杂的理论分析, 但是其结果具有随机性, 精度很难提高: 为了增加一位小数点的精度,即误差减小到原来的 \(\frac{1}{10}\) ,重复试验次数需要增加到原来的 100 倍。 随机模拟方法有如下特点:

{\color{red}\textbf{[量级估计]:} 这就是为什么蒙特卡洛方法``计算密集''。想要从2位精度($N=10^4$)提升到4位精度,需要$N=10^8$——计算量增加一万倍!但好处是这个收敛速度与问题维度无关,高维积分时反而比传统数值方法更有优势。}

\begin{itemize}
\item 应用面广,适应性强。只要问题能够清楚地描述出来, 就可以编写模拟程序产生大量数据, 通过分析模拟数据解决问题。
\end{itemize}

\begin{itemize}
\item 算法简单, 容易改变条件, 但计算量大。
\end{itemize}

\begin{itemize}
\item 结果具有随机性,精度较低 (一般为 \(O\left( \frac{1}{\sqrt{N}}\right)\) 级)。
\end{itemize}

\begin{itemize}
\item 模拟结果的收敛服从概率论规律。
\end{itemize}

\begin{itemize}
\item 对维数增加不敏感。在计算定积分时, 如果使用传统数值算法, 维数增加会造成计算时间指数增加,但是如果使用随机模拟方法计算,则维数增加仅仅造成不多的影响。
\end{itemize}

{\color{red}\textbf{[为什么对维度不敏感？]:} 传统数值积分像``网格扫描''——在每个维度上取$n$个点,总点数是$n^d$(维度灾难)。蒙特卡洛像``随机投点''——无论多少维,取$N$个随机点,误差都是$O(1/\sqrt{N})$。这是蒙特卡洛在高维问题(如量子多体、金融衍生品定价)中不可替代的原因。}

随机模拟用在科学研究中, 常常作为探索性试验来使用。假设科学家有了一个新的模型或技术的想法, 但是不知道它的效果怎样所以还没有对其进行深入的理论分析, 就可以用随机模拟方法大量地重复生成模拟数据, 根据多次重复的总体效果来判断这种模型或技术的性能。如果模拟获得了好的结果, 再进行深入理论分析对模型进行完善; 如果模拟发现了这个模型的缺点, 可以进行有针对性的修改, 或者考虑转而其它解决办法。

{\color{red}\textbf{[研究策略]:} 蒙特卡洛模拟是``理论想法的试验场''。就像工程师先做纸板模型测试,再投入生产——花几小时写模拟代码,跑几天计算,就能知道新想法是``金矿''还是``死胡同'',避免在错误方向上浪费数月理论推导时间。}

随机模拟在科学研究中的另一种作用是说明新的模型或技术的有效性。在公开发表的统计学论文中, 已经有一半以上的文章包括随机模拟结果 (也叫数值结果), 用来辅助说明自己提出的模型或方法的有效性。有时因为对模型或方法很难进行彻底的理论分析, 仅仅使用大量的随机模拟结果来说明模型或方法的有效性。当然, 因为模型都是有可变参数的, 随机模拟只能针对某些参数组合给出结果, 所以, 一般认为仅有模拟结果而没有理论分析结果的研究论文是不全面的。§3.5给出一个用随机模拟说明统计技术优良性的例子。

除了以上应用, 随机模拟还是许多新的统计方法的主要工具, 例如, 蒙特卡洛检验, bootstrap 置信区间和 bootstrap 偏差修正, MCMC。利用大量计算机计算 (包括随机模拟) 来进行统计推断的统计学分支叫做 ``计算统计'' (computational statistics), 在本章后面各节将介绍随机模拟的一些应用和技巧。

{\color{red}\textbf{[展望]:} 从1940年代诞生到今天,蒙特卡洛方法已从核武器设计的保密工具,演变为现代科学计算的基础设施——从统计推断(MCMC采样贝叶斯后验)到机器学习(dropout训练神经网络),从量子化学到金融工程,无处不在。计算能力每提升100倍,蒙特卡洛方法的应用边界就扩展10倍。}
